Statistical analysis of plasma carotenoids, plasma cytokines and immune cells from randomized cross-over USDA inflammation clinical trial. Subjects consumed both low lycopene tomato (yellow) and high lycopene tomato-soy juices (red) for 4 weeks each.
Crossover clinical trial design supplementing individuals with obesity 360 mL of a low carotenoid tomato juice or a high lycopene tomato-soy juice daily. Daily serving of low carotenoid tomato juice consisted of ~1.5mg lycopene/day while high lycopene tomato-soy juice intervention consisted of 54 mg lycopene/day in addition to 210 mg total soy isoflavones/day.
library(tidyverse) # data wrangling
library(readxl) # read in excel files
library(janitor) # clean up names in dataset
library(corrr) # finding correlations
library(rstatix) # stats
library(lme4) # mixed linear modeling
library(knitr) # aesthetic table viewing
library(lmerTest) # add pvalue column to lmer models
library(purrr) # create functions
library(broom.mixed) # generate tidy data frames for lmer results
library(MuMIn) # lmer model testing using AICc
library(kableExtra)
library(ggthemes)
library(ggtext)
library(ggpubr)
# load data
meta_table <- read_excel("CompiledData_Results_Meta.xlsx",
sheet = "metadata_corrected_withsequence")
# clean up variable names
meta_table <- clean_names(meta_table)
str(meta_table)
## tibble [60 × 89] (S3: tbl_df/tbl/data.frame)
## $ patient_id : num [1:60] 6101 6102 6103 6104 6105 ...
## $ sex : chr [1:60] "F" "M" "F" "M" ...
## $ age_at_enrollment : num [1:60] 58 65 60 54 40 57 68 75 36 46 ...
## $ bmi_at_enrollment : num [1:60] 31.1 36.8 30.3 33.1 30.4 ...
## $ intervention : chr [1:60] "Red" "Red" "Red" "Yellow" ...
## $ sequence : chr [1:60] "R_Y" "R_Y" "R_Y" "Y_R" ...
## $ intervention_week : num [1:60] 2 2 2 2 2 2 2 2 2 2 ...
## $ pre_post : chr [1:60] "pre" "pre" "pre" "pre" ...
## $ period : chr [1:60] "B1" "B1" "B1" "B1" ...
## $ if_ng : num [1:60] 1.16 11.74 0.93 0.63 1.88 ...
## $ il_1b : num [1:60] 15.26 28.75 8.62 11.04 22.5 ...
## $ il_2 : chr [1:60] "0.98" "9.32" "0.55000000000000004" "0.65" ...
## $ il_6 : num [1:60] 2.18 5.83 2.53 1.29 2.25 ...
## $ il_8 : num [1:60] 2.86 4.35 3.93 1.93 4.14 3.05 1.35 2.01 4.43 2.86 ...
## $ il_10 : chr [1:60] "0.05" "0.71" "1.89" "7.17" ...
## $ il_12p70 : num [1:60] 2.86 60.93 5.46 2.79 4.06 ...
## $ mcp_1 : num [1:60] 51.9 189.3 224.9 244.7 167.9 ...
## $ tn_fa : num [1:60] 14.9 92.3 26.5 19.6 49.6 ...
## $ il_13 : chr [1:60] "9.64" "92.04" "21.61" "23.17" ...
## $ il_5 : num [1:60] 3 13.04 2.74 5.45 5.38 ...
## $ il_1ra : num [1:60] 3.52 11.36 3.71 3.38 6.08 ...
## $ il_12p40 : num [1:60] 51.5 139.4 81 78.7 175.9 ...
## $ gm_csf : num [1:60] 10.9 195.7 17.9 27.4 72.4 ...
## $ il_4 : chr [1:60] "NA" "8.8000000000000007" "1.6" "4.49" ...
## $ x01_cd45_cd66b_lymph_dc_mono : num [1:60] 99.3 95.6 99.7 99.7 96.4 ...
## $ x02_cd45_cd66b_grans : num [1:60] 0.682 0.33 0.305 0.285 0.696 ...
## $ cd45_cd14 : num [1:60] 46.9 71 69.8 55.8 42.4 ...
## $ cd45_cd20 : num [1:60] 44.2 59.9 60.6 54 31.9 ...
## $ x03_cd3_cd45_cd3_t_cells : num [1:60] 25.7 43.7 54.5 14.2 14.3 ...
## $ x04_tc_rgd_cd3_ab_t_cells : num [1:60] 22 42.6 53.1 11.2 13.6 ...
## $ x05_cd4_cd8_cd8_t_cells : num [1:60] 11.96 3.84 13.23 6.5 5.72 ...
## $ ccr7_cd8 : num [1:60] 4.385 1.54 4.186 0.962 1.882 ...
## $ ccr7_cd8_2 : num [1:60] 7.56 2.28 9.04 5.52 3.83 ...
## $ x06_cd45ro_cd45ra_naive_cd8 : num [1:60] 3.692 0.233 3.753 0.652 1.667 ...
## $ x07_cd46ro_cd45ra_cm_cd8 : num [1:60] 0.506 0.955 0.263 0.196 0.146 ...
## $ x08_cd45ro_cd45ra_em_cd8 : num [1:60] 6.555 0.395 7.505 4.528 1.807 ...
## $ x09_cd45r0_cd45ra_te_cd8 : num [1:60] 0.826 1.322 1.2 0.775 1.661 ...
## $ x10_cd38_hladr_activated_cd8 : num [1:60] 0.2543 0.0823 0.1523 0.1347 0.0742 ...
## $ x11_cd4_cd8_cd4_t_cells : num [1:60] 9.14 37.12 37.51 3.64 6.98 ...
## $ ccr7_cd4 : num [1:60] 7.98 33.92 32.87 2.49 5.3 ...
## $ ccr7_cd4_2 : num [1:60] 1.16 3.19 4.64 1.15 1.68 ...
## $ x12_cd45ro_cd45ra_naive_cd4 : num [1:60] 3.163 3.77 17.23 0.693 1.731 ...
## $ x13_cd45ro_cd45ra_cm_cd4 : num [1:60] 3.51 22.39 8.5 1.25 2.63 ...
## $ x14_cd45ro_cd45ra : num [1:60] 1.053 2.964 3.919 0.893 1.488 ...
## $ x15_cd45ro_cd45ra_te_cd4 : num [1:60] 0.0257 0.0664 0.5457 0.1154 0.092 ...
## $ x16_cd38_hladr_activated_cd4 : num [1:60] 0.321 0.698 0.384 0.123 0.14 ...
## $ cd45_cxcr5_th : num [1:60] 8.291 28.987 0.133 2.95 6.456 ...
## $ x17_cd25_cd127_tregs : num [1:60] 0.348 0.545 0.438 0.131 0.193 ...
## $ x18_ccr4_cd4_total_ccr4_treg : num [1:60] 0.317 0.494 0.316 0.123 0.18 ...
## $ x19_cd45ra_cd45ro_ccr4_treg_naive : num [1:60] 0.00973 0.03494 0.00676 0.00955 0.00744 ...
## $ x21_cd45ra_cd45ro_ccr4_treg_memory : num [1:60] 0.3017 0.4366 0.2905 0.0949 0.166 ...
## $ x20_hladr_total_ccr4_treg_activated : num [1:60] 0.2867 0.2547 0.1017 0.0838 0.1119 ...
## $ x22_cxcr3_ccr6_th1 : num [1:60] 0.7981 2.0038 0.0209 0.4805 0.8362 ...
## $ x23_cxcr3_ccr6_th2 : num [1:60] 5.6202 13.4081 0.0747 1.2236 4.3881 ...
## $ x24_cxcr3_ccr6_th17 : num [1:60] 1.1118 9.5315 0.0153 0.6605 0.3389 ...
## $ cd16_cd161 : num [1:60] 27.1 46.9 55.1 15 23.5 ...
## $ x25_cd19_cd3_b_cells : num [1:60] 2.53 8.22 7.86 1.62 9.48 ...
## $ cd20 : num [1:60] 0.115 0.111 0.2 0.147 0.102 ...
## $ cd20_2 : num [1:60] 2.41 8.09 7.64 1.47 9.38 ...
## $ x26_cd27_ig_d_naive_b_cells : num [1:60] 1.19 6.99 4.62 1.29 6.69 ...
## $ x27_cd27_ig_d_memory_b_cells : num [1:60] 0.697 0.529 2.22 0.101 0.878 ...
## $ x28_cd27_ig_d_memory_resting_b_cells: num [1:60] 0.1083 0.2849 0.304 0.0499 1.4241 ...
## $ x30_cd27_cd38_plasmablasts : num [1:60] 0.0634 0.0277 0.1582 0.0378 0.062 ...
## $ cd19_cd3 : num [1:60] 66.8 35.6 30.2 79.9 65.8 ...
## $ x31_cd14_monocytes : num [1:60] 47.7 19.5 23.2 41.8 48.7 ...
## $ x32_cd16_non_classical_mono : num [1:60] 6.24 3.29 3.07 7.15 3.4 ...
## $ x33_cd16_classical_mono : num [1:60] 38.5 16.1 19.1 33.5 45.1 ...
## $ x34_hladr_cd56 : num [1:60] 43.8 16.7 20.3 33.9 41.1 ...
## $ cd16_cd123 : num [1:60] 1.142 0.272 0.657 0.645 0.613 ...
## $ x35_cd16_cd123_cd11c_p_dc : num [1:60] 0.709 0.248 0.594 0.37 0.495 ...
## $ cd16_cd123_2 : num [1:60] 33.1 10.5 15.1 25.6 37.3 ...
## $ x36_cd16_cd123_cd11c_m_dc : num [1:60] 27.7 10.4 14.7 23.6 37.1 ...
## $ cd123 : num [1:60] 65.7 35.1 29.5 73.5 65.1 ...
## $ total_dc : num [1:60] 29.5 10.9 15.9 34 38.2 ...
## $ x37_cd56_cd161_cd123_nk_cells : num [1:60] 13.01 9.43 2.44 5.33 14.75 ...
## $ x38_cd16_nk_cells : num [1:60] 1.601 1.552 0.874 28.656 5.754 ...
## $ cd16_nk_cells : num [1:60] 11.41 7.85 1.56 41.75 8.9 ...
## $ cd14_cd11b : num [1:60] 46.4 17.3 20.5 35 48.7 ...
## $ x40_cd14_mdsc_mono : num [1:60] 28.983 7.399 9.73 0.215 39.928 ...
## $ cd66b_cd11b : num [1:60] 0.314 0.251 0.265 0.195 0.532 ...
## $ x41_cd66b_mdsc_grans : num [1:60] 0.2057 0.0724 0.0985 NA 0.4165 ...
## $ lutein : num [1:60] 735 120 556 419 240 ...
## $ zeaxanthin : num [1:60] 94.6 140.1 99 81.6 53.7 ...
## $ b_cryptoxanthin : num [1:60] 314 286 389 230 570 ...
## $ a_carotene : num [1:60] 195.7 101.8 181.2 112.6 37.1 ...
## $ b_carotene : num [1:60] 927 309 434 536 121 ...
## $ other_cis_lyc : num [1:60] 151 276 161 165 273 ...
## $ all_trans_lyc : num [1:60] 115 259 185 175 336 ...
## $ x5_cis_lyc : num [1:60] 132 262 166 149 235 ...
# convert variables that should be factors to factors
meta_table <- meta_table %>%
mutate(across(.cols = c("patient_id", "period",
"intervention", "intervention_week",
"pre_post", "sex", "sequence"),
.fns = as.factor))
# some stuff came in as characters but should be numeric
meta_table <- meta_table %>%
mutate(across(.cols = c("il_2", "il_10", "il_13", "il_4"),
.fns = as.numeric))
str(meta_table)
## tibble [60 × 89] (S3: tbl_df/tbl/data.frame)
## $ patient_id : Factor w/ 12 levels "6101","6102",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : Factor w/ 2 levels "F","M": 1 2 1 2 2 2 2 1 2 2 ...
## $ age_at_enrollment : num [1:60] 58 65 60 54 40 57 68 75 36 46 ...
## $ bmi_at_enrollment : num [1:60] 31.1 36.8 30.3 33.1 30.4 ...
## $ intervention : Factor w/ 3 levels "Baseline","Red",..: 2 2 2 3 2 3 2 3 3 3 ...
## $ sequence : Factor w/ 2 levels "R_Y","Y_R": 1 1 1 2 1 2 1 2 2 2 ...
## $ intervention_week : Factor w/ 5 levels "0","2","6","10",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ pre_post : Factor w/ 3 levels "baseline","post",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ period : Factor w/ 5 levels "B0","B1","B2",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ if_ng : num [1:60] 1.16 11.74 0.93 0.63 1.88 ...
## $ il_1b : num [1:60] 15.26 28.75 8.62 11.04 22.5 ...
## $ il_2 : num [1:60] 0.98 9.32 0.55 0.65 1.8 0.35 0.38 0.05 NA 0.16 ...
## $ il_6 : num [1:60] 2.18 5.83 2.53 1.29 2.25 ...
## $ il_8 : num [1:60] 2.86 4.35 3.93 1.93 4.14 3.05 1.35 2.01 4.43 2.86 ...
## $ il_10 : num [1:60] 0.05 0.71 1.89 7.17 0.75 2.2 1.91 4.26 0.02 0.48 ...
## $ il_12p70 : num [1:60] 2.86 60.93 5.46 2.79 4.06 ...
## $ mcp_1 : num [1:60] 51.9 189.3 224.9 244.7 167.9 ...
## $ tn_fa : num [1:60] 14.9 92.3 26.5 19.6 49.6 ...
## $ il_13 : num [1:60] 9.64 92.04 21.61 23.17 60.22 ...
## $ il_5 : num [1:60] 3 13.04 2.74 5.45 5.38 ...
## $ il_1ra : num [1:60] 3.52 11.36 3.71 3.38 6.08 ...
## $ il_12p40 : num [1:60] 51.5 139.4 81 78.7 175.9 ...
## $ gm_csf : num [1:60] 10.9 195.7 17.9 27.4 72.4 ...
## $ il_4 : num [1:60] NA 8.8 1.6 4.49 1.22 1.12 0.31 0.91 0.07 NA ...
## $ x01_cd45_cd66b_lymph_dc_mono : num [1:60] 99.3 95.6 99.7 99.7 96.4 ...
## $ x02_cd45_cd66b_grans : num [1:60] 0.682 0.33 0.305 0.285 0.696 ...
## $ cd45_cd14 : num [1:60] 46.9 71 69.8 55.8 42.4 ...
## $ cd45_cd20 : num [1:60] 44.2 59.9 60.6 54 31.9 ...
## $ x03_cd3_cd45_cd3_t_cells : num [1:60] 25.7 43.7 54.5 14.2 14.3 ...
## $ x04_tc_rgd_cd3_ab_t_cells : num [1:60] 22 42.6 53.1 11.2 13.6 ...
## $ x05_cd4_cd8_cd8_t_cells : num [1:60] 11.96 3.84 13.23 6.5 5.72 ...
## $ ccr7_cd8 : num [1:60] 4.385 1.54 4.186 0.962 1.882 ...
## $ ccr7_cd8_2 : num [1:60] 7.56 2.28 9.04 5.52 3.83 ...
## $ x06_cd45ro_cd45ra_naive_cd8 : num [1:60] 3.692 0.233 3.753 0.652 1.667 ...
## $ x07_cd46ro_cd45ra_cm_cd8 : num [1:60] 0.506 0.955 0.263 0.196 0.146 ...
## $ x08_cd45ro_cd45ra_em_cd8 : num [1:60] 6.555 0.395 7.505 4.528 1.807 ...
## $ x09_cd45r0_cd45ra_te_cd8 : num [1:60] 0.826 1.322 1.2 0.775 1.661 ...
## $ x10_cd38_hladr_activated_cd8 : num [1:60] 0.2543 0.0823 0.1523 0.1347 0.0742 ...
## $ x11_cd4_cd8_cd4_t_cells : num [1:60] 9.14 37.12 37.51 3.64 6.98 ...
## $ ccr7_cd4 : num [1:60] 7.98 33.92 32.87 2.49 5.3 ...
## $ ccr7_cd4_2 : num [1:60] 1.16 3.19 4.64 1.15 1.68 ...
## $ x12_cd45ro_cd45ra_naive_cd4 : num [1:60] 3.163 3.77 17.23 0.693 1.731 ...
## $ x13_cd45ro_cd45ra_cm_cd4 : num [1:60] 3.51 22.39 8.5 1.25 2.63 ...
## $ x14_cd45ro_cd45ra : num [1:60] 1.053 2.964 3.919 0.893 1.488 ...
## $ x15_cd45ro_cd45ra_te_cd4 : num [1:60] 0.0257 0.0664 0.5457 0.1154 0.092 ...
## $ x16_cd38_hladr_activated_cd4 : num [1:60] 0.321 0.698 0.384 0.123 0.14 ...
## $ cd45_cxcr5_th : num [1:60] 8.291 28.987 0.133 2.95 6.456 ...
## $ x17_cd25_cd127_tregs : num [1:60] 0.348 0.545 0.438 0.131 0.193 ...
## $ x18_ccr4_cd4_total_ccr4_treg : num [1:60] 0.317 0.494 0.316 0.123 0.18 ...
## $ x19_cd45ra_cd45ro_ccr4_treg_naive : num [1:60] 0.00973 0.03494 0.00676 0.00955 0.00744 ...
## $ x21_cd45ra_cd45ro_ccr4_treg_memory : num [1:60] 0.3017 0.4366 0.2905 0.0949 0.166 ...
## $ x20_hladr_total_ccr4_treg_activated : num [1:60] 0.2867 0.2547 0.1017 0.0838 0.1119 ...
## $ x22_cxcr3_ccr6_th1 : num [1:60] 0.7981 2.0038 0.0209 0.4805 0.8362 ...
## $ x23_cxcr3_ccr6_th2 : num [1:60] 5.6202 13.4081 0.0747 1.2236 4.3881 ...
## $ x24_cxcr3_ccr6_th17 : num [1:60] 1.1118 9.5315 0.0153 0.6605 0.3389 ...
## $ cd16_cd161 : num [1:60] 27.1 46.9 55.1 15 23.5 ...
## $ x25_cd19_cd3_b_cells : num [1:60] 2.53 8.22 7.86 1.62 9.48 ...
## $ cd20 : num [1:60] 0.115 0.111 0.2 0.147 0.102 ...
## $ cd20_2 : num [1:60] 2.41 8.09 7.64 1.47 9.38 ...
## $ x26_cd27_ig_d_naive_b_cells : num [1:60] 1.19 6.99 4.62 1.29 6.69 ...
## $ x27_cd27_ig_d_memory_b_cells : num [1:60] 0.697 0.529 2.22 0.101 0.878 ...
## $ x28_cd27_ig_d_memory_resting_b_cells: num [1:60] 0.1083 0.2849 0.304 0.0499 1.4241 ...
## $ x30_cd27_cd38_plasmablasts : num [1:60] 0.0634 0.0277 0.1582 0.0378 0.062 ...
## $ cd19_cd3 : num [1:60] 66.8 35.6 30.2 79.9 65.8 ...
## $ x31_cd14_monocytes : num [1:60] 47.7 19.5 23.2 41.8 48.7 ...
## $ x32_cd16_non_classical_mono : num [1:60] 6.24 3.29 3.07 7.15 3.4 ...
## $ x33_cd16_classical_mono : num [1:60] 38.5 16.1 19.1 33.5 45.1 ...
## $ x34_hladr_cd56 : num [1:60] 43.8 16.7 20.3 33.9 41.1 ...
## $ cd16_cd123 : num [1:60] 1.142 0.272 0.657 0.645 0.613 ...
## $ x35_cd16_cd123_cd11c_p_dc : num [1:60] 0.709 0.248 0.594 0.37 0.495 ...
## $ cd16_cd123_2 : num [1:60] 33.1 10.5 15.1 25.6 37.3 ...
## $ x36_cd16_cd123_cd11c_m_dc : num [1:60] 27.7 10.4 14.7 23.6 37.1 ...
## $ cd123 : num [1:60] 65.7 35.1 29.5 73.5 65.1 ...
## $ total_dc : num [1:60] 29.5 10.9 15.9 34 38.2 ...
## $ x37_cd56_cd161_cd123_nk_cells : num [1:60] 13.01 9.43 2.44 5.33 14.75 ...
## $ x38_cd16_nk_cells : num [1:60] 1.601 1.552 0.874 28.656 5.754 ...
## $ cd16_nk_cells : num [1:60] 11.41 7.85 1.56 41.75 8.9 ...
## $ cd14_cd11b : num [1:60] 46.4 17.3 20.5 35 48.7 ...
## $ x40_cd14_mdsc_mono : num [1:60] 28.983 7.399 9.73 0.215 39.928 ...
## $ cd66b_cd11b : num [1:60] 0.314 0.251 0.265 0.195 0.532 ...
## $ x41_cd66b_mdsc_grans : num [1:60] 0.2057 0.0724 0.0985 NA 0.4165 ...
## $ lutein : num [1:60] 735 120 556 419 240 ...
## $ zeaxanthin : num [1:60] 94.6 140.1 99 81.6 53.7 ...
## $ b_cryptoxanthin : num [1:60] 314 286 389 230 570 ...
## $ a_carotene : num [1:60] 195.7 101.8 181.2 112.6 37.1 ...
## $ b_carotene : num [1:60] 927 309 434 536 121 ...
## $ other_cis_lyc : num [1:60] 151 276 161 165 273 ...
## $ all_trans_lyc : num [1:60] 115 259 185 175 336 ...
## $ x5_cis_lyc : num [1:60] 132 262 166 149 235 ...
# changing factor levels for pre_post
meta_table$pre_post <- factor(meta_table$pre_post,
levels = c("pre", "post"))
levels(meta_table$pre_post)
## [1] "pre" "post"
# Calculate total_cis_lyc, total_lyc, and total_carotenoids
meta_table <- meta_table %>%
rename(n5_cis_lyc = x5_cis_lyc) %>%
mutate(total_cis_lyc = other_cis_lyc + n5_cis_lyc,
total_lyc = all_trans_lyc + total_cis_lyc,
total_carotenoids = lutein + zeaxanthin + b_cryptoxanthin +
a_carotene + b_carotene + total_lyc)
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = all_trans_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_bw() +
labs(x = "Intervention Week",
y = "All-trans-lycopene levels (nmol/L)",
title = "All-trans-lycopene levels in each patient before/after each intervention")
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = total_cis_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id)) +
theme_bw() +
labs(x = "Intervention Week",
y = "Total cis-lycopene levels (nmol/L)",
title = "Total cis-lycopene levels in each patient before/after each intervention")
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = total_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_bw() +
labs(x = "Intervention Week",
y = "Total lycopene levels (nmol/L)",
title = "Total lycopene levels in each patient before/after each intervention")
# create a more specific pre_post_intervention column
meta_table_edited <- meta_table %>%
unite(col = "pre_post_intervention",
c("pre_post","intervention"),
sep = "_",
remove = FALSE)
# make pre_post_intervention column factors
meta_table_edited$pre_post_intervention <- as.factor(meta_table_edited$pre_post_intervention)
# relevel factor columns
meta_table_edited$pre_post_intervention <- factor(meta_table_edited$pre_post_intervention, levels = c("pre_Yellow", "post_Yellow", "pre_Red", "post_Red"))
meta_table_edited$intervention <- factor(meta_table_edited$intervention,
levels = c("Yellow", "Red"))
# make legend title
legendtitle_ppintervention <- "Timepoint"
# labels
labs_ppintervention <- c("pre control",
"post control",
"pre Tomato-Soy",
"post Tomato-Soy")
meta_table_edited %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = intervention, y = total_lyc, fill = pre_post_intervention)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(legendtitle_ppintervention,
values = c("pre_Red" = "#FF9966",
"post_Red" = "#FF3300",
"pre_Yellow" = "#FFFF99",
"post_Yellow" = "yellow"),
labels = labs_ppintervention) +
theme_clean() +
labs(x = "",
y = "Total lycopene levels (nmol/L)",
title = "Total lycopene levels before/after juice interventions")
meta_table_edited %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = pre_post_intervention, y = total_lyc, fill = intervention)) +
geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1")) +
scale_x_discrete(labels = labs_ppintervention) +
theme_clean(base_size = 22, base_family = "sans") +
labs(x = "",
y = "Overall plasma lycopene conc. (nmol/L)",
title = "Total plasma lycopene levels",
subtitle = "Before and after juice interventions")
# publish-ready plot
library(ggpubr)
(total_lyc_levels <- meta_table_edited %>%
filter(intervention != "Baseline") %>%
ggpaired(x = "pre_post", y = "total_lyc", fill = "intervention", line.color = "gray", line.size = 1, facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", linewidth = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "nmol/L plasma",
title = "Concentration of Lycopene",
subtitle = ""))
# convert meta_table_edited to long format for lycopene
meta_table_lyc_long <- meta_table_edited %>%
pivot_longer(cols = total_lyc,
names_to = "total_lycopene",
values_to = "nmol_per_L")
lyc_long_subset <- meta_table_lyc_long %>%
select(patient_id, nmol_per_L, pre_post_intervention) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = nmol_per_L) %>%
mutate(red_FC = post_Red/pre_Red,
yellow_FC = post_Yellow/pre_Yellow)
lyc_long_subset %>%
summarize(mean_red_FC = mean(red_FC),
mean_yellow_FC = mean(yellow_FC))
## # A tibble: 1 × 2
## mean_red_FC mean_yellow_FC
## <dbl> <dbl>
## 1 2.83 0.966
Plot histogram
gghistogram(meta_table_lyc_long$nmol_per_L, bins = 40)
Shapiro’s normality test
# shapiro normality test for total lycopene
meta_table_lyc_long %>%
group_by(intervention) %>%
shapiro_test(vars = "nmol_per_L")
## # A tibble: 3 × 4
## intervention variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Yellow nmol_per_L 0.876 0.00698
## 2 Red nmol_per_L 0.848 0.00202
## 3 <NA> nmol_per_L 0.913 0.233
P val for shapiro test turned out to < 0.05 for both interventions, suggesting data here is not normal. Let’s log transform before running paired t-tests.
meta_table_lyc_log2_long <- meta_table_lyc_long %>%
mutate(log2_nmol_per_L = log2(nmol_per_L))
Without transforming, here’s a paired t-test
compare_means(nmol_per_L ~ pre_post, meta_table_lyc_log2_long, method = "t.test", paired = TRUE, group.by = "intervention")
## # A tibble: 2 × 9
## intervention .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Red nmol_per_L pre post 0.00104 0.0021 0.001 ** T-test
## 2 Yellow nmol_per_L pre post 0.709 0.71 0.709 ns T-test
Plot histogram
gghistogram(meta_table_lyc_log2_long$log2_nmol_per_L, bins = 40)
Shapiro’s normality test
# shapiro normality test for total lycopene
meta_table_lyc_log2_long %>%
group_by(intervention) %>%
shapiro_test(vars = "log2_nmol_per_L")
## # A tibble: 3 × 4
## intervention variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Yellow log2_nmol_per_L 0.928 0.0869
## 2 Red log2_nmol_per_L 0.964 0.521
## 3 <NA> log2_nmol_per_L 0.903 0.173
compare_means(log2_nmol_per_L ~ pre_post, meta_table_lyc_log2_long, method = "t.test", paired = TRUE, group.by = "intervention")
## # A tibble: 2 × 9
## intervention .y. group1 group2 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Red log2_nmol… pre post 3.84e-5 7.70e-5 3.8e-05 **** T-test
## 2 Yellow log2_nmol… pre post 3.21e-1 3.2 e-1 0.32 ns T-test
Total lycopene levels are significantly increasing only after post-Red intervention
# convert cytokines from wide to long
cytokines_long <- meta_table[-c(25:81)] %>%
pivot_longer(cols = if_ng:il_4,
names_to = "cytokines",
values_to = "cyto_conc_pg_ml")
#cytokines_long_nona <- na.omit(cytokines_long)
note 11/30/23: I took this part out since NAs don’t matter for lmm
After testing several models, the best model turned out to have pre_post as a fixed variable and patient_id as random effect. Carotenoids added no statistically significant effect to the model, suggesting that carotenoids are not contributing to cytokine levels. Set REML = FALSE when comparing models.
# turn off sci notation
options(scipen = 999)
# sequence effect test lmer model function
cyto_seq_function <- function(cytokines_long) {
lmer(cyto_conc_pg_ml ~ sequence + (1|patient_id),
data = cytokines_long,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_seq_models <- map(split(cytokines_long,
cytokines_long$cytokines),
cyto_seq_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_seq_results <- map_df(cyto_seq_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_seq_results_coefonly <- cyto_seq_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
# results
kable(cyto_seq_results_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | sequenceY_R | -27.444 | 24.681 | -1.112 | 10.000 | 0.292 | ns |
| if_ng | fixed | NA | sequenceY_R | -1.194 | 1.477 | -0.808 | 10.000 | 0.438 | ns |
| il_10 | fixed | NA | sequenceY_R | 1.728 | 1.220 | 1.416 | 9.996 | 0.187 | ns |
| il_12p40 | fixed | NA | sequenceY_R | -16.115 | 23.825 | -0.676 | 10.000 | 0.514 | ns |
| il_12p70 | fixed | NA | sequenceY_R | -7.289 | 7.609 | -0.958 | 10.000 | 0.361 | ns |
| il_13 | fixed | NA | sequenceY_R | -18.773 | 13.629 | -1.377 | 10.724 | 0.196 | ns |
| il_1b | fixed | NA | sequenceY_R | -1.552 | 8.358 | -0.186 | 10.000 | 0.856 | ns |
| il_1ra | fixed | NA | sequenceY_R | -2.149 | 1.685 | -1.275 | 10.000 | 0.231 | ns |
| il_2 | fixed | NA | sequenceY_R | -1.585 | 1.440 | -1.101 | 8.191 | 0.302 | ns |
| il_4 | fixed | NA | sequenceY_R | -0.667 | 1.778 | -0.375 | 6.000 | 0.721 | ns |
| il_5 | fixed | NA | sequenceY_R | 1.657 | 3.027 | 0.547 | 10.000 | 0.596 | ns |
| il_6 | fixed | NA | sequenceY_R | -2.432 | 1.896 | -1.283 | 10.000 | 0.229 | ns |
| il_8 | fixed | NA | sequenceY_R | 0.432 | 0.674 | 0.641 | 10.000 | 0.536 | ns |
| mcp_1 | fixed | NA | sequenceY_R | 49.779 | 28.416 | 1.752 | 10.000 | 0.110 | ns |
| tn_fa | fixed | NA | sequenceY_R | -6.730 | 10.112 | -0.666 | 10.000 | 0.521 | ns |
Are there any significant cytokines?
# extract statistically significant cytokines
cytokines_seq_psig <- cyto_seq_results_coefonly %>%
filter(cyto_seq_results_coefonly$p.value < 0.05)
print(cytokines_seq_psig$cytokines)
## character(0)
# filter data for post-intervention only
cytokines_post_long <- cytokines_long %>%
filter(intervention != "Baseline") %>%
filter(pre_post == "post")
# treatment effect test lmer model function
cyto_trt_function <- function(cytokines_post_long) {
lmer(cyto_conc_pg_ml ~ intervention + (1|patient_id),
data = cytokines_post_long,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trt_models <- map(split(cytokines_post_long,
cytokines_post_long$cytokines),
cyto_trt_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trt_results <- map_df(cyto_trt_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_trt_coefonly <- cyto_trt_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
kable(cyto_trt_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | interventionYellow | -2.237 | 2.904 | -0.770 | 11.000 | 0.457 | ns |
| if_ng | fixed | NA | interventionYellow | 0.081 | 0.213 | 0.380 | 11.000 | 0.711 | ns |
| il_10 | fixed | NA | interventionYellow | 0.121 | 0.158 | 0.766 | 10.035 | 0.461 | ns |
| il_12p40 | fixed | NA | interventionYellow | 6.302 | 5.595 | 1.126 | 11.000 | 0.284 | ns |
| il_12p70 | fixed | NA | interventionYellow | -0.220 | 0.361 | -0.609 | 11.000 | 0.555 | ns |
| il_13 | fixed | NA | interventionYellow | -0.213 | 1.716 | -0.124 | 7.124 | 0.905 | ns |
| il_1b | fixed | NA | interventionYellow | 1.167 | 1.119 | 1.043 | 11.000 | 0.319 | ns |
| il_1ra | fixed | NA | interventionYellow | 0.311 | 0.469 | 0.662 | 11.000 | 0.521 | ns |
| il_2 | fixed | NA | interventionYellow | -0.153 | 0.168 | -0.909 | 6.131 | 0.398 | ns |
| il_4 | fixed | NA | interventionYellow | -0.129 | 0.075 | -1.713 | 7.000 | 0.131 | ns |
| il_5 | fixed | NA | interventionYellow | 0.500 | 0.708 | 0.706 | 11.000 | 0.495 | ns |
| il_6 | fixed | NA | interventionYellow | 0.991 | 0.750 | 1.322 | 11.000 | 0.213 | ns |
| il_8 | fixed | NA | interventionYellow | -0.060 | 0.355 | -0.169 | 11.000 | 0.869 | ns |
| mcp_1 | fixed | NA | interventionYellow | 3.200 | 8.404 | 0.381 | 11.000 | 0.711 | ns |
| tn_fa | fixed | NA | interventionYellow | 0.458 | 1.097 | 0.418 | 11.000 | 0.684 | ns |
# extract statistically significant cytokines
cytokines_trt_psig <- cyto_trt_coefonly %>%
filter(cyto_trt_coefonly$p.value < 0.05)
print(cytokines_trt_psig$cytokines)
## character(0)
# subset data into yellow results
cytokines_Y <- cytokines_long %>%
filter(intervention == "Yellow")
# treatment effect test lmer model function
cyto_trtY_function <- function(cytokines_Y) {
lmer(cyto_conc_pg_ml ~ pre_post + (1|patient_id),
data = cytokines_Y,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtY_models <- map(split(cytokines_Y,
cytokines_Y$cytokines),
cyto_trtY_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtY_results <- map_df(cyto_trtY_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_trtY_coefonly <- cyto_trtY_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
kable(cyto_trtY_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | pre_postpost | -10.773 | 9.940 | -1.084 | 11.000 | 0.302 | ns |
| if_ng | fixed | NA | pre_postpost | -0.582 | 0.554 | -1.050 | 11.000 | 0.316 | ns |
| il_10 | fixed | NA | pre_postpost | 0.003 | 0.144 | 0.024 | 9.060 | 0.982 | ns |
| il_12p40 | fixed | NA | pre_postpost | -7.124 | 9.540 | -0.747 | 11.000 | 0.471 | ns |
| il_12p70 | fixed | NA | pre_postpost | -3.552 | 3.510 | -1.012 | 11.000 | 0.333 | ns |
| il_13 | fixed | NA | pre_postpost | -10.171 | 8.663 | -1.174 | 7.959 | 0.274 | ns |
| il_1b | fixed | NA | pre_postpost | -4.913 | 3.071 | -1.600 | 11.000 | 0.138 | ns |
| il_1ra | fixed | NA | pre_postpost | -1.743 | 1.465 | -1.189 | 11.000 | 0.259 | ns |
| il_2 | fixed | NA | pre_postpost | -1.097 | 0.952 | -1.151 | 7.491 | 0.285 | ns |
| il_4 | fixed | NA | pre_postpost | -0.752 | 0.722 | -1.043 | 7.000 | 0.332 | ns |
| il_5 | fixed | NA | pre_postpost | -0.906 | 1.413 | -0.641 | 11.000 | 0.535 | ns |
| il_6 | fixed | NA | pre_postpost | -0.614 | 0.659 | -0.932 | 11.000 | 0.372 | ns |
| il_8 | fixed | NA | pre_postpost | 0.130 | 0.296 | 0.439 | 11.000 | 0.669 | ns |
| mcp_1 | fixed | NA | pre_postpost | 1.432 | 5.523 | 0.259 | 11.000 | 0.800 | ns |
| tn_fa | fixed | NA | pre_postpost | -3.615 | 4.206 | -0.859 | 11.000 | 0.408 | ns |
Significant cytokines
# extract statistically significant cytokines
cytokines_trtY_psig <- cyto_trtY_coefonly %>%
filter(cyto_trtY_coefonly$p.value < 0.05)
print(cytokines_trtY_psig$cytokines)
## character(0)
# subset data into yellow results
cytokines_R <- cytokines_long %>%
filter(intervention == "Red")
# treatment effect test lmer model function
cyto_trtR_function <- function(cytokines_R) {
lmer(cyto_conc_pg_ml ~ pre_post + (1|patient_id),
data = cytokines_R,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtR_models <- map(split(cytokines_R,
cytokines_R$cytokines),
cyto_trtR_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtR_results <- map_df(cyto_trtR_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_trtR_coefonly <- cyto_trtR_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
kable(cyto_trtR_coefonly, format = "markdown", digits = 3) %>%
row_spec(11, bold = TRUE)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | pre_postpost | -7.572 | 5.759 | -1.315 | 11.000 | 0.215 | ns |
| if_ng | fixed | NA | pre_postpost | -0.497 | 0.412 | -1.209 | 11.000 | 0.252 | ns |
| il_10 | fixed | NA | pre_postpost | -0.157 | 0.139 | -1.131 | 10.028 | 0.284 | ns |
| il_12p40 | fixed | NA | pre_postpost | -13.373 | 6.652 | -2.010 | 11.000 | 0.070 | ns |
| il_12p70 | fixed | NA | pre_postpost | -3.119 | 2.415 | -1.291 | 11.000 | 0.223 | ns |
| il_13 | fixed | NA | pre_postpost | -5.090 | 4.748 | -1.072 | 7.669 | 0.316 | ns |
| il_1b | fixed | NA | pre_postpost | -1.908 | 1.306 | -1.461 | 11.000 | 0.172 | ns |
| il_1ra | fixed | NA | pre_postpost | -0.435 | 0.453 | -0.961 | 11.000 | 0.357 | ns |
| il_2 | fixed | NA | pre_postpost | -0.522 | 0.556 | -0.938 | 7.313 | 0.378 | ns |
| il_4 | fixed | NA | pre_postpost | -0.581 | 0.527 | -1.104 | 7.000 | 0.306 | ns |
| il_5 | fixed | NA | pre_postpost | -1.808 | 0.815 | -2.220 | 11.000 | 0.048 | * |
| il_6 | fixed | NA | pre_postpost | -1.102 | 0.820 | -1.343 | 11.000 | 0.206 | ns |
| il_8 | fixed | NA | pre_postpost | 0.073 | 0.217 | 0.338 | 11.000 | 0.741 | ns |
| mcp_1 | fixed | NA | pre_postpost | -2.799 | 6.896 | -0.406 | 11.000 | 0.693 | ns |
| tn_fa | fixed | NA | pre_postpost | -4.860 | 2.934 | -1.656 | 11.000 | 0.126 | ns |
Significant cytokines
# extract statistically significant cytokines
cytokines_trtR_psig <- cyto_trtR_coefonly %>%
filter(cyto_trtR_coefonly$p.value < 0.05)
print(cytokines_trtR_psig$cytokines)
## [1] "il_5"
meta_table %>%
filter(intervention == "Red") %>%
ggplot(aes(x = intervention_week, y = il_5,
color = intervention, group = intervention)) +
geom_line() +
theme_classic() +
labs(x = "Week",
y = "IL-5 levels (pg/mL)",
title = "IL-5 changes by subject before and after red intervention",
color = "Patient ID") +
facet_wrap(~patient_id,
scales = "free_y") +
scale_color_manual(values = c("Red" = "tomato1"))
meta_table %>%
filter(intervention == "Red") %>%
ggplot(aes(x = pre_post, y = il_5, color = patient_id)) +
geom_line(aes(group = patient_id)) +
labs(x = "",
y = "IL-5 conc (pg/mL)",
title = "IL-5 levels in each patient before and after red intervention") +
theme_classic()
meta_table %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = pre_post, y = il_5, color = intervention)) +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_classic() +
labs(x = "",
y = "IL-5 conc (pg/mL)",
title = "IL-5 levels in each patient pre- and post- red and yellow intervention")
cytokines_long %>%
filter(intervention == "Red") %>%
filter(cytokines == "il_5") %>%
ggplot(aes(x = pre_post, y = cyto_conc_pg_ml, fill = intervention)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter() +
scale_fill_manual(values = c("Red" = "tomato1")) +
labs(x = "",
y = "IL- 5 levels (pg/mL)",
title = "IL-5 levels before and after red intervention") +
theme_bw()
I know the data isn’t normally distributed. Mixed linear models do not assume normality of the dependant variable, however the residuals should be normally distributed.
I need to figure out how to make a function that extracts makes a QQplot for each model. But since IL-5 was significant, let’s look at how the residuals look for this cytokine (for every comparison) before log transforming
#qqplot_fx <- function(cytokine){
#qqnorm(resid(cyto_seq_models$cytokine))
#qqline(resid(cyto_seq_models$cytokine))
#}
#qqplot_fx("il_5")
#cyto_seq_results$
qqnorm(resid(cyto_seq_models$il_5))
qqline(resid(cyto_seq_models$il_5))
That doesn’t look good.
qqnorm(resid(cyto_trt_models$il_5))
qqline(resid(cyto_trt_models$il_5))
doesn’t look as bad as sequence
qqnorm(resid(cyto_trtR_models$il_5))
qqline(resid(cyto_trtR_models$il_5))
doesn’t look that bad for Red
qqnorm(resid(cyto_trtY_models$il_5))
qqline(resid(cyto_trtY_models$il_5))
I will log transform the data and perform mixed linear modeling and paired t-tests for every comparison to see how the data compares.
Wrangle
cytokines_log2conc_long <- cytokines_long %>%
mutate(log2_cyto_conc_pg_ml = log2(cyto_conc_pg_ml))
We want to make sure there are no seq effects
# sequence effect test lmer model function
cyto_log2_seq_function <- function(cytokines_log2conc_long) {
lmer(log2_cyto_conc_pg_ml ~ sequence + (1|patient_id),
data = cytokines_log2conc_long,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_log2_seq_models <- map(split(cytokines_log2conc_long,
cytokines_log2conc_long$cytokines),
cyto_log2_seq_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_log2_seq_results <- map_df(cyto_log2_seq_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_log2_seq_results_coefonly <- cyto_log2_seq_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
# before log transform
qqnorm(resid(cyto_seq_models$il_5))
qqline(resid(cyto_seq_models$il_5))
# after log transform
qqnorm(resid(cyto_log2_seq_models$il_5))
qqline(resid(cyto_log2_seq_models$il_5))
This looks better after log transform.
# results
kable(cyto_log2_seq_results_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | sequenceY_R | -0.703 | 0.813 | -0.864 | 10.000 | 0.408 | ns |
| if_ng | fixed | NA | sequenceY_R | -0.639 | 1.065 | -0.600 | 10.000 | 0.562 | ns |
| il_10 | fixed | NA | sequenceY_R | 0.897 | 1.380 | 0.650 | 9.855 | 0.530 | ns |
| il_12p40 | fixed | NA | sequenceY_R | -0.181 | 0.408 | -0.444 | 10.000 | 0.666 | ns |
| il_12p70 | fixed | NA | sequenceY_R | -0.804 | 0.889 | -0.905 | 10.000 | 0.387 | ns |
| il_13 | fixed | NA | sequenceY_R | -1.081 | 0.848 | -1.275 | 10.230 | 0.230 | ns |
| il_1b | fixed | NA | sequenceY_R | -0.321 | 0.836 | -0.384 | 10.000 | 0.709 | ns |
| il_1ra | fixed | NA | sequenceY_R | -0.570 | 0.518 | -1.101 | 10.000 | 0.297 | ns |
| il_2 | fixed | NA | sequenceY_R | -1.546 | 1.196 | -1.292 | 8.035 | 0.232 | ns |
| il_4 | fixed | NA | sequenceY_R | -0.557 | 1.491 | -0.374 | 6.000 | 0.721 | ns |
| il_5 | fixed | NA | sequenceY_R | 0.115 | 0.595 | 0.194 | 10.000 | 0.850 | ns |
| il_6 | fixed | NA | sequenceY_R | -0.976 | 0.607 | -1.607 | 10.000 | 0.139 | ns |
| il_8 | fixed | NA | sequenceY_R | 0.287 | 0.371 | 0.773 | 10.000 | 0.457 | ns |
| mcp_1 | fixed | NA | sequenceY_R | 0.469 | 0.304 | 1.543 | 10.000 | 0.154 | ns |
| tn_fa | fixed | NA | sequenceY_R | -0.033 | 0.404 | -0.083 | 10.000 | 0.936 | ns |
Are there any significant cytokines?
# extract statistically significant cytokines
cytokines_log2_seq_psig <- cyto_log2_seq_results_coefonly %>%
filter(cyto_log2_seq_results_coefonly$p.value < 0.05)
print(cytokines_log2_seq_psig$cytokines)
## character(0)
# filter data for post-intervention only
cytokines_log2conc_post_long <- cytokines_log2conc_long %>%
filter(intervention != "Baseline") %>%
filter(pre_post == "post")
# treatment effect test lmer model function
cyto_log2_trt_function <- function(cytokines_log2conc_post_long) {
lmer(log2_cyto_conc_pg_ml ~ intervention + (1|patient_id),
data = cytokines_log2conc_post_long,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_log2_trt_models <- map(split(cytokines_log2conc_post_long,
cytokines_log2conc_post_long$cytokines),
cyto_log2_trt_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_log2_trt_results <- map_df(cyto_log2_trt_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_log2_trt_coefonly <- cyto_log2_trt_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
# before log transform
qqnorm(resid(cyto_trt_models$il_5))
qqline(resid(cyto_trt_models$il_5))
# log transformed
qqnorm(resid(cyto_log2_trt_models$il_5))
qqline(resid(cyto_log2_trt_models$il_5))
Don’t know if this helps for intervention comparison
kable(cyto_log2_trt_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | interventionYellow | 0.116 | 0.118 | 0.980 | 11.000 | 0.348 | ns |
| if_ng | fixed | NA | interventionYellow | 0.161 | 0.230 | 0.700 | 11.000 | 0.498 | ns |
| il_10 | fixed | NA | interventionYellow | 0.183 | 0.110 | 1.665 | 9.985 | 0.127 | ns |
| il_12p40 | fixed | NA | interventionYellow | 0.177 | 0.153 | 1.159 | 11.000 | 0.271 | ns |
| il_12p70 | fixed | NA | interventionYellow | 0.166 | 0.153 | 1.083 | 11.000 | 0.302 | ns |
| il_13 | fixed | NA | interventionYellow | 0.146 | 0.146 | 0.997 | 7.041 | 0.352 | ns |
| il_1b | fixed | NA | interventionYellow | 0.126 | 0.164 | 0.767 | 11.000 | 0.459 | ns |
| il_1ra | fixed | NA | interventionYellow | 0.083 | 0.155 | 0.535 | 11.000 | 0.603 | ns |
| il_2 | fixed | NA | interventionYellow | 0.216 | 0.291 | 0.745 | 6.079 | 0.484 | ns |
| il_4 | fixed | NA | interventionYellow | -0.052 | 0.070 | -0.743 | 7.000 | 0.482 | ns |
| il_5 | fixed | NA | interventionYellow | 0.379 | 0.216 | 1.749 | 11.000 | 0.108 | ns |
| il_6 | fixed | NA | interventionYellow | 0.331 | 0.263 | 1.260 | 11.000 | 0.234 | ns |
| il_8 | fixed | NA | interventionYellow | -0.105 | 0.156 | -0.674 | 11.000 | 0.514 | ns |
| mcp_1 | fixed | NA | interventionYellow | 0.007 | 0.066 | 0.108 | 11.000 | 0.916 | ns |
| tn_fa | fixed | NA | interventionYellow | 0.035 | 0.071 | 0.499 | 11.000 | 0.627 | ns |
# extract statistically significant cytokines
cytokines_log2_trt_psig <- cyto_log2_trt_coefonly %>%
filter(cyto_log2_trt_coefonly$p.value < 0.05)
print(cytokines_log2_trt_psig$cytokines)
## character(0)
# subset data into yellow results
cytokines_Y_log2 <- cytokines_log2conc_long %>%
filter(intervention == "Yellow")
# treatment effect test lmer model function
cyto_trtY_log2_function <- function(cytokines_Y_log2) {
lmer(log2_cyto_conc_pg_ml ~ pre_post + (1|patient_id),
data = cytokines_Y_log2,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtY_log2_models <- map(split(cytokines_Y_log2,
cytokines_Y_log2$cytokines),
cyto_trtY_log2_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtY_log2_results <- map_df(cyto_trtY_log2_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_trtY_log2_coefonly <- cyto_trtY_log2_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
# before log transform
qqnorm(resid(cyto_trtY_models$il_5))
qqline(resid(cyto_trtY_models$il_5))
# after log transform
qqnorm(resid(cyto_trtY_log2_models$il_5))
qqline(resid(cyto_trtY_log2_models$il_5))
kable(cyto_trtY_log2_coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | pre_postpost | -0.183 | 0.135 | -1.353 | 11.000 | 0.203 | ns |
| if_ng | fixed | NA | pre_postpost | -0.235 | 0.165 | -1.423 | 11.000 | 0.182 | ns |
| il_10 | fixed | NA | pre_postpost | -0.241 | 0.226 | -1.066 | 8.920 | 0.314 | ns |
| il_12p40 | fixed | NA | pre_postpost | -0.068 | 0.126 | -0.543 | 11.000 | 0.598 | ns |
| il_12p70 | fixed | NA | pre_postpost | -0.053 | 0.145 | -0.365 | 11.000 | 0.722 | ns |
| il_13 | fixed | NA | pre_postpost | -0.130 | 0.242 | -0.539 | 6.933 | 0.607 | ns |
| il_1b | fixed | NA | pre_postpost | -0.206 | 0.208 | -0.991 | 11.000 | 0.343 | ns |
| il_1ra | fixed | NA | pre_postpost | -0.213 | 0.220 | -0.970 | 11.000 | 0.353 | ns |
| il_2 | fixed | NA | pre_postpost | -0.029 | 0.338 | -0.087 | 7.055 | 0.933 | ns |
| il_4 | fixed | NA | pre_postpost | -0.121 | 0.192 | -0.633 | 7.000 | 0.547 | ns |
| il_5 | fixed | NA | pre_postpost | 0.124 | 0.239 | 0.518 | 11.000 | 0.615 | ns |
| il_6 | fixed | NA | pre_postpost | -0.125 | 0.356 | -0.350 | 11.000 | 0.733 | ns |
| il_8 | fixed | NA | pre_postpost | -0.024 | 0.141 | -0.170 | 11.000 | 0.868 | ns |
| mcp_1 | fixed | NA | pre_postpost | -0.012 | 0.043 | -0.278 | 11.000 | 0.786 | ns |
| tn_fa | fixed | NA | pre_postpost | -0.082 | 0.094 | -0.876 | 11.000 | 0.400 | ns |
Significant cytokines
# extract statistically significant cytokines
cytokines_trtY_log2_psig <- cyto_trtY_log2_coefonly %>%
filter(cyto_trtY_log2_coefonly$p.value < 0.05)
print(cytokines_trtY_log2_psig$cytokines)
## character(0)
# subset data into yellow results
cytokines_R_log2 <- cytokines_log2conc_long %>%
filter(intervention == "Red")
# treatment effect test lmer model function
cyto_trtR_log2_function <- function(cytokines_R_log2) {
lmer(log2_cyto_conc_pg_ml ~ pre_post + (1|patient_id),
data = cytokines_R_log2,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cyto_trtR_log2_models <- map(split(cytokines_R_log2,
cytokines_R_log2$cytokines),
cyto_trtR_log2_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cyto_trtR_log2_results <- map_df(cyto_trtR_log2_models,
tidy,
.id = "cytokines")
# extract fixed coefficients for sequence only
cyto_trtR_log2coefonly <- cyto_trtR_log2_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
# before log transform
qqnorm(resid(cyto_trtR_models$il_5))
qqline(resid(cyto_trtR_models$il_5))
# log transformed
qqnorm(resid(cyto_trtR_log2_models$il_5))
qqline(resid(cyto_trtR_log2_models$il_5))
kable(cyto_trtR_log2coefonly, format = "markdown", digits = 3)
| cytokines | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| gm_csf | fixed | NA | pre_postpost | -0.308 | 0.103 | -2.985 | 11.000 | 0.012 | * |
| if_ng | fixed | NA | pre_postpost | -0.307 | 0.197 | -1.559 | 11.000 | 0.147 | ns |
| il_10 | fixed | NA | pre_postpost | -0.137 | 0.422 | -0.325 | 9.926 | 0.752 | ns |
| il_12p40 | fixed | NA | pre_postpost | -0.288 | 0.156 | -1.840 | 11.000 | 0.093 | ns |
| il_12p70 | fixed | NA | pre_postpost | -0.622 | 0.240 | -2.596 | 11.000 | 0.025 | * |
| il_13 | fixed | NA | pre_postpost | -0.130 | 0.199 | -0.652 | 7.181 | 0.535 | ns |
| il_1b | fixed | NA | pre_postpost | -0.227 | 0.164 | -1.387 | 11.000 | 0.193 | ns |
| il_1ra | fixed | NA | pre_postpost | -0.109 | 0.120 | -0.912 | 11.000 | 0.381 | ns |
| il_2 | fixed | NA | pre_postpost | -0.273 | 0.277 | -0.988 | 7.142 | 0.356 | ns |
| il_4 | fixed | NA | pre_postpost | -0.200 | 0.156 | -1.282 | 7.000 | 0.241 | ns |
| il_5 | fixed | NA | pre_postpost | -0.373 | 0.172 | -2.166 | 11.000 | 0.053 | ns |
| il_6 | fixed | NA | pre_postpost | -0.200 | 0.161 | -1.240 | 11.000 | 0.241 | ns |
| il_8 | fixed | NA | pre_postpost | 0.068 | 0.080 | 0.842 | 11.000 | 0.418 | ns |
| mcp_1 | fixed | NA | pre_postpost | 0.006 | 0.047 | 0.126 | 11.000 | 0.902 | ns |
| tn_fa | fixed | NA | pre_postpost | -0.174 | 0.084 | -2.070 | 11.000 | 0.063 | ns |
Significant cytokines
# extract statistically significant cytokines
cytokines_trtR_log2_psig <- cyto_trtR_log2coefonly %>%
filter(cyto_trtR_log2coefonly$p.value < 0.05)
print(cytokines_trtR_log2_psig$cytokines)
## [1] "gm_csf" "il_12p70"
note 11/30/23: Log transforming changes results a bit. IL-5 is no longer significant (pval = 0.053), and now GM-CSF and IL-12p70 are significant. I’ve done nonparametric stats on untransformed data, these 2 cytokines along with IL-5 were significant. The question is, should I be doing log transformation for the mixed linear models? Q-Q plot looks better for IL-5 after log transforming…
compare_means(log2_cyto_conc_pg_ml ~ pre_post, cytokines_R_log2, method = "t.test", paired = TRUE, group.by = "cytokines")
## # A tibble: 15 × 9
## cytokines .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 if_ng log2_cyto_conc… pre post 0.147 1 0.147 ns T-test
## 2 il_1b log2_cyto_conc… pre post 0.193 1 0.193 ns T-test
## 3 il_2 log2_cyto_conc… pre post 0.350 1 0.350 ns T-test
## 4 il_6 log2_cyto_conc… pre post 0.241 1 0.241 ns T-test
## 5 il_8 log2_cyto_conc… pre post 0.418 1 0.418 ns T-test
## 6 il_10 log2_cyto_conc… pre post 0.847 1 0.847 ns T-test
## 7 il_12p70 log2_cyto_conc… pre post 0.0249 0.35 0.025 * T-test
## 8 mcp_1 log2_cyto_conc… pre post 0.902 1 0.902 ns T-test
## 9 tn_fa log2_cyto_conc… pre post 0.0628 0.75 0.063 ns T-test
## 10 il_13 log2_cyto_conc… pre post 0.540 1 0.540 ns T-test
## 11 il_5 log2_cyto_conc… pre post 0.0531 0.69 0.053 ns T-test
## 12 il_1ra log2_cyto_conc… pre post 0.381 1 0.381 ns T-test
## 13 il_12p40 log2_cyto_conc… pre post 0.0928 1 0.093 ns T-test
## 14 gm_csf log2_cyto_conc… pre post 0.0124 0.19 0.012 * T-test
## 15 il_4 log2_cyto_conc… pre post 0.241 1 0.241 ns T-test
We get the same result for paired t-test on log transformed data
note 11/30/23: if there are any NAs, should I assume below LOD? Or were values not recorded?
# convert immune cell data from wide to long
cells_long <- meta_table %>%
select(1:24, starts_with("x")) %>%
pivot_longer(cols = x01_cd45_cd66b_lymph_dc_mono:x41_cd66b_mdsc_grans,
names_to = "cell_type",
values_to = "cell_value")
# lmer model function
cells_seq_function <- function(cells_long) {
lmer(cell_value ~ sequence + (1|patient_id), data = cells_long, REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cells_seq_models <- map(split(cells_long,
cells_long$cell_type),
cells_seq_function)
# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_seq_results <- map_df(cells_seq_models,
tidy,
.id = "cell_type")
cells_seq_results_coef <- cells_seq_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
kable(cells_seq_results_coef, format = "markdown", digits = 3)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | sequenceY_R | -0.024 | 2.152 | -0.011 | 10.000 | 0.991 | ns |
| x02_cd45_cd66b_grans | fixed | NA | sequenceY_R | 0.445 | 0.558 | 0.797 | 10.000 | 0.444 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | sequenceY_R | 5.038 | 7.441 | 0.677 | 10.000 | 0.514 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | sequenceY_R | 4.702 | 7.682 | 0.612 | 10.000 | 0.554 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | sequenceY_R | 3.751 | 2.477 | 1.514 | 10.000 | 0.161 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | sequenceY_R | 1.338 | 1.832 | 0.730 | 10.000 | 0.482 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | sequenceY_R | 0.152 | 0.210 | 0.724 | 10.000 | 0.485 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | sequenceY_R | 1.094 | 1.512 | 0.723 | 10.000 | 0.486 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | sequenceY_R | -0.109 | 0.437 | -0.249 | 10.000 | 0.809 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | sequenceY_R | -0.002 | 0.026 | -0.089 | 10.000 | 0.931 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | sequenceY_R | 0.546 | 7.572 | 0.072 | 10.000 | 0.944 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | sequenceY_R | 1.412 | 4.877 | 0.290 | 10.000 | 0.778 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | sequenceY_R | -0.649 | 3.423 | -0.190 | 10.000 | 0.853 | ns |
| x14_cd45ro_cd45ra | fixed | NA | sequenceY_R | 0.835 | 0.737 | 1.133 | 10.000 | 0.284 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | sequenceY_R | -0.337 | 0.565 | -0.597 | 10.001 | 0.564 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | sequenceY_R | 0.001 | 0.087 | 0.011 | 10.000 | 0.991 | ns |
| x17_cd25_cd127_tregs | fixed | NA | sequenceY_R | 0.290 | 0.296 | 0.978 | 10.000 | 0.351 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | sequenceY_R | 0.312 | 0.268 | 1.163 | 10.000 | 0.272 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | sequenceY_R | 0.019 | 0.011 | 1.802 | 10.000 | 0.102 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | sequenceY_R | 0.008 | 0.061 | 0.137 | 10.000 | 0.893 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | sequenceY_R | 0.224 | 0.204 | 1.097 | 10.000 | 0.298 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | sequenceY_R | 0.634 | 0.708 | 0.896 | 10.000 | 0.391 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | sequenceY_R | 7.379 | 5.644 | 1.307 | 10.000 | 0.220 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | sequenceY_R | 0.684 | 2.080 | 0.329 | 10.000 | 0.749 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | sequenceY_R | 1.080 | 2.731 | 0.395 | 10.000 | 0.701 | ns |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | sequenceY_R | 1.249 | 2.409 | 0.518 | 10.000 | 0.615 | ns |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | sequenceY_R | -0.236 | 0.322 | -0.733 | 10.000 | 0.480 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | sequenceY_R | -0.177 | 0.219 | -0.806 | 10.000 | 0.439 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | sequenceY_R | -0.033 | 0.021 | -1.613 | 10.000 | 0.138 | ns |
| x31_cd14_monocytes | fixed | NA | sequenceY_R | -5.481 | 6.350 | -0.863 | 10.000 | 0.408 | ns |
| x32_cd16_non_classical_mono | fixed | NA | sequenceY_R | -0.440 | 0.703 | -0.627 | 10.000 | 0.545 | ns |
| x33_cd16_classical_mono | fixed | NA | sequenceY_R | -4.477 | 5.667 | -0.790 | 10.000 | 0.448 | ns |
| x34_hladr_cd56 | fixed | NA | sequenceY_R | -4.413 | 5.412 | -0.815 | 10.000 | 0.434 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | sequenceY_R | -0.092 | 0.069 | -1.346 | 10.000 | 0.208 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | sequenceY_R | -3.230 | 4.565 | -0.708 | 10.000 | 0.495 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | sequenceY_R | -1.323 | 2.273 | -0.582 | 10.000 | 0.573 | ns |
| x38_cd16_nk_cells | fixed | NA | sequenceY_R | 1.153 | 2.027 | 0.569 | 10.000 | 0.582 | ns |
| x40_cd14_mdsc_mono | fixed | NA | sequenceY_R | -4.353 | 3.288 | -1.324 | 10.000 | 0.215 | ns |
| x41_cd66b_mdsc_grans | fixed | NA | sequenceY_R | -0.035 | 0.128 | -0.274 | 10.000 | 0.789 | ns |
# extract statistically significant cytokines
cells_seq_psig <- cells_seq_results_coef %>%
filter(cells_seq_results_coef$p.value < 0.05)
print(cells_seq_psig$cell_type)
## character(0)
*No sequence effects
# filter data for post-intervention only
cells_post_long <- cells_long %>%
filter(intervention != "Baseline") %>%
filter(pre_post == "post")
str(cells_post_long)
## tibble [936 × 26] (S3: tbl_df/tbl/data.frame)
## $ patient_id : Factor w/ 12 levels "6101","6102",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
## $ age_at_enrollment: num [1:936] 58 58 58 58 58 58 58 58 58 58 ...
## $ bmi_at_enrollment: num [1:936] 31.1 31.1 31.1 31.1 31.1 ...
## $ intervention : Factor w/ 3 levels "Baseline","Red",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ sequence : Factor w/ 2 levels "R_Y","Y_R": 1 1 1 1 1 1 1 1 1 1 ...
## $ intervention_week: Factor w/ 5 levels "0","2","6","10",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ pre_post : Factor w/ 2 levels "pre","post": 2 2 2 2 2 2 2 2 2 2 ...
## $ period : Factor w/ 5 levels "B0","B1","B2",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ if_ng : num [1:936] 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 ...
## $ il_1b : num [1:936] 17.9 17.9 17.9 17.9 17.9 ...
## $ il_2 : num [1:936] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ...
## $ il_6 : num [1:936] 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 ...
## $ il_8 : num [1:936] 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 3.53 ...
## $ il_10 : num [1:936] 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 ...
## $ il_12p70 : num [1:936] 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 2.79 ...
## $ mcp_1 : num [1:936] 60.4 60.4 60.4 60.4 60.4 ...
## $ tn_fa : num [1:936] 16.8 16.8 16.8 16.8 16.8 ...
## $ il_13 : num [1:936] 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 7.47 ...
## $ il_5 : num [1:936] 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 3.74 ...
## $ il_1ra : num [1:936] 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 4.29 ...
## $ il_12p40 : num [1:936] 45.3 45.3 45.3 45.3 45.3 ...
## $ gm_csf : num [1:936] 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 6.46 ...
## $ il_4 : num [1:936] NA NA NA NA NA NA NA NA NA NA ...
## $ cell_type : chr [1:936] "x01_cd45_cd66b_lymph_dc_mono" "x02_cd45_cd66b_grans" "x03_cd3_cd45_cd3_t_cells" "x04_tc_rgd_cd3_ab_t_cells" ...
## $ cell_value : num [1:936] 98.61 1.36 36.48 31.92 15.05 ...
# remove x41_cd66b_mdsc_grans because it's causing an error
cells_post_long_edited <- cells_post_long %>%
filter(cell_type != "x41_cd66b_mdsc_grans")
# treatment effect test lmer model function
cells_trt_function <- function(cells_post_long_edited) {
lmer(cell_value ~ intervention + (1|patient_id),
data = cells_post_long_edited,
REML = TRUE)
}
# break data up into subsets based on cell, then apply funtion to each subset
cells_trt_models <- map(split(cells_post_long_edited,
cells_post_long_edited$cell_type),
cells_trt_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trt_results <- map_df(cells_trt_models,
tidy,
.id = "cell_type")
# extract fixed coefficients for sequence only
cells_trt_coefonly <- cells_trt_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
kable(cells_trt_coefonly, format = "markdown", digits = 4)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | interventionYellow | 0.1437 | 1.5446 | 0.0930 | 11 | 0.9276 | ns |
| x02_cd45_cd66b_grans | fixed | NA | interventionYellow | -0.3867 | 0.3053 | -1.2669 | 11 | 0.2313 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | interventionYellow | -1.9095 | 4.0782 | -0.4682 | 11 | 0.6488 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | interventionYellow | -2.0179 | 4.0621 | -0.4968 | 11 | 0.6291 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | interventionYellow | -0.0913 | 1.2534 | -0.0729 | 11 | 0.9432 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | interventionYellow | 0.1903 | 0.6862 | 0.2773 | 11 | 0.7867 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | interventionYellow | -0.0137 | 0.0717 | -0.1914 | 11 | 0.8517 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | interventionYellow | -0.6942 | 1.1263 | -0.6164 | 11 | 0.5502 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | interventionYellow | 0.3789 | 0.2313 | 1.6379 | 11 | 0.1297 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | interventionYellow | 0.0168 | 0.0119 | 1.4113 | 11 | 0.1858 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | interventionYellow | -1.6574 | 2.9437 | -0.5630 | 11 | 0.5847 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | interventionYellow | -0.6033 | 1.9806 | -0.3046 | 11 | 0.7664 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | interventionYellow | -0.1943 | 0.9081 | -0.2140 | 11 | 0.8345 | ns |
| x14_cd45ro_cd45ra | fixed | NA | interventionYellow | 0.0149 | 0.2516 | 0.0594 | 11 | 0.9537 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | interventionYellow | -1.0865 | 1.1069 | -0.9816 | 22 | 0.3370 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | interventionYellow | 0.0301 | 0.0376 | 0.8014 | 11 | 0.4398 | ns |
| x17_cd25_cd127_tregs | fixed | NA | interventionYellow | 0.1232 | 0.0949 | 1.2986 | 11 | 0.2207 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | interventionYellow | 0.1147 | 0.0958 | 1.1970 | 11 | 0.2565 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | interventionYellow | 0.0058 | 0.0069 | 0.8350 | 11 | 0.4215 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | interventionYellow | 0.0527 | 0.0390 | 1.3536 | 11 | 0.2030 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | interventionYellow | 0.0754 | 0.0581 | 1.2975 | 11 | 0.2210 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | interventionYellow | -0.5504 | 0.6386 | -0.8618 | 11 | 0.4072 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | interventionYellow | -0.4578 | 1.7276 | -0.2650 | 11 | 0.7959 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | interventionYellow | 1.2526 | 1.1705 | 1.0701 | 11 | 0.3075 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | interventionYellow | 0.1449 | 0.8671 | 0.1671 | 11 | 0.8703 | ns |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | interventionYellow | 0.0560 | 0.8292 | 0.0676 | 11 | 0.9473 | ns |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | interventionYellow | 0.0936 | 0.2156 | 0.4343 | 11 | 0.6725 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | interventionYellow | -0.0552 | 0.0831 | -0.6645 | 11 | 0.5200 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | interventionYellow | 0.0185 | 0.0134 | 1.3743 | 11 | 0.1967 | ns |
| x31_cd14_monocytes | fixed | NA | interventionYellow | 0.4091 | 3.8171 | 0.1072 | 11 | 0.9166 | ns |
| x32_cd16_non_classical_mono | fixed | NA | interventionYellow | -0.3330 | 0.3102 | -1.0735 | 11 | 0.3060 | ns |
| x33_cd16_classical_mono | fixed | NA | interventionYellow | 0.6321 | 3.5425 | 0.1784 | 11 | 0.8616 | ns |
| x34_hladr_cd56 | fixed | NA | interventionYellow | 0.6282 | 3.2256 | 0.1947 | 11 | 0.8491 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | interventionYellow | 0.0082 | 0.0508 | 0.1606 | 11 | 0.8753 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | interventionYellow | 0.6105 | 2.9571 | 0.2064 | 11 | 0.8402 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | interventionYellow | -4.8544 | 1.6409 | -2.9583 | 11 | 0.0130 | * |
| x38_cd16_nk_cells | fixed | NA | interventionYellow | 3.3508 | 0.8126 | 4.1235 | 11 | 0.0017 | ** |
| x40_cd14_mdsc_mono | fixed | NA | interventionYellow | -10.9242 | 3.2788 | -3.3317 | 11 | 0.0067 | ** |
# extract statistically significant cell type
cells_trt_psig <- cells_trt_coefonly %>%
filter(cells_trt_coefonly$p.value < 0.05)
print(cells_trt_psig$cell_type)
## [1] "x37_cd56_cd161_cd123_nk_cells" "x38_cd16_nk_cells"
## [3] "x40_cd14_mdsc_mono"
*3 significant cell types
sigcells_trt <- cells_trt_psig$cell_type
cells_post_long_edited %>%
filter(cell_type %in% sigcells_trt) %>%
ggplot(aes(x = intervention, y = cell_value, color = patient_id)) +
geom_line(aes(group = patient_id)) +
labs(x = "",
y = "",
title = "") +
theme_classic() +
facet_wrap(~ cell_type, scales = "free_y")
cells_post_long_edited %>%
filter(cell_type %in% sigcells_trt) %>%
ggplot(aes(x = intervention, y = cell_value, fill = intervention)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
geom_jitter() +
labs(x = "Intervention",
y = "Cell value",
title = "Significantly different cell types: post- yellow vs. post- red interventions") +
facet_wrap(~ cell_type, ncol = 3, nrow = 1, scales = "free_y") +
theme_bw()
x38_CD16 NK cells are significantly higher in post- yellow compared to post-red
x37_cd56_cd161_cd123_nk_cells and x40_cd14_mdsc_mono significantly lower post-yellow intervention compared to post-red
# subset data into yellow results
cells_Y_long <- cells_long %>%
filter(intervention == "Yellow")
# remove x41_cd66b_mdsc_grans because it has too many NAs
cells_Y_long_edited <- cells_Y_long %>%
filter(cell_type != "x41_cd66b_mdsc_grans")
# treatment effect test lmer model function
cells_trtY_function <- function(cells_Y_long_edited) {
lmer(cell_value ~ pre_post + (1|patient_id),
data = cells_Y_long_edited,
REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cells_trtY_models <- map(split(cells_Y_long_edited,
cells_Y_long_edited$cell_type),
cells_trtY_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trtY_results <- map_df(cells_trtY_models,
tidy,
.id = "cell_type")
# extract fixed coefficients for sequence only
cells_trtY_coefonly <- cells_trtY_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
# extract statistically significant cytokines
cells_trtY_psig <- cells_trtY_coefonly %>%
filter(cells_trtY_coefonly$p.value < 0.05)
print(cells_trtY_psig$cell_type)
## [1] "x05_cd4_cd8_cd8_t_cells" "x08_cd45ro_cd45ra_em_cd8"
## [3] "x09_cd45r0_cd45ra_te_cd8" "x15_cd45ro_cd45ra_te_cd4"
## [5] "x26_cd27_ig_d_naive_b_cells"
sigcells_Y <- cells_trtY_psig$cell_type
cells_Y_long_edited %>%
filter(cell_type %in% sigcells_Y) %>%
ggplot(aes(x = pre_post, y = cell_value, color = patient_id)) +
geom_line(aes(group = patient_id)) +
labs(x = "",
y = "",
title = "") +
theme_classic() +
facet_wrap(~ cell_type, scales = "free_y")
cells_Y_long_edited %>%
filter(cell_type %in% sigcells_Y) %>%
ggplot(aes(x = pre_post, y = cell_value, fill = pre_post)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(values = c("pre" = "gray",
"post" = "gold")) +
geom_jitter() +
labs(x = "Pre vs post",
y = "Cell value",
title = "Significantly different cell types: pre vs post yellow") +
facet_wrap(~ cell_type, ncol = 3, nrow = 2, scales = "free_y") +
theme_bw()
*5 significant cell types
# subset data into red results
cells_R_long <- cells_long %>%
filter(intervention == "Red")
# remove x41_cd66b_mdsc_grans because it has too many NAs
cells_R_long_edited <- cells_R_long %>%
filter(cell_type != "x41_cd66b_mdsc_grans")
# treatment effect test lmer model function
cells_trtR_function <- function(cells_R_long_edited) {
lmer(cell_value ~ pre_post + (1|patient_id),
data = cells_R_long_edited,
REML = TRUE)
}
# break data up into subsets based on cell, then apply funtion to each subset
cells_trtR_models <- map(split(cells_R_long_edited,
cells_R_long_edited$cell_type),
cells_trtR_function)
# create single data frame of results by applying the lmer function to each element in list. apply tidy from broom.mixed package to clean up results
cells_trtR_results <- map_df(cells_trtR_models,
tidy,
.id = "cell_type")
# extract fixed coefficients for sequence only
cells_trtR_coefonly <- cells_trtR_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
Results
kable(cells_trtR_coefonly, format = "markdown", digits = 3)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | pre_postpost | -0.476 | 0.978 | -0.487 | 11 | 0.636 | ns |
| x02_cd45_cd66b_grans | fixed | NA | pre_postpost | 0.354 | 0.228 | 1.556 | 11 | 0.148 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | pre_postpost | 0.294 | 6.161 | 0.048 | 11 | 0.963 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | pre_postpost | 0.086 | 6.144 | 0.014 | 11 | 0.989 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | pre_postpost | 1.983 | 1.198 | 1.655 | 11 | 0.126 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | pre_postpost | 0.137 | 0.463 | 0.296 | 11 | 0.773 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | pre_postpost | 0.063 | 0.049 | 1.296 | 11 | 0.221 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | pre_postpost | 1.645 | 1.099 | 1.496 | 11 | 0.163 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | pre_postpost | 0.147 | 0.191 | 0.770 | 11 | 0.457 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | pre_postpost | -0.005 | 0.014 | -0.346 | 11 | 0.736 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | pre_postpost | -2.279 | 5.045 | -0.452 | 11 | 0.660 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | pre_postpost | -1.932 | 2.696 | -0.717 | 11 | 0.488 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | pre_postpost | -0.585 | 0.924 | -0.633 | 11 | 0.540 | ns |
| x14_cd45ro_cd45ra | fixed | NA | pre_postpost | 0.374 | 0.255 | 1.467 | 11 | 0.170 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | pre_postpost | 1.147 | 1.105 | 1.037 | 22 | 0.311 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | pre_postpost | -0.041 | 0.031 | -1.304 | 11 | 0.219 | ns |
| x17_cd25_cd127_tregs | fixed | NA | pre_postpost | -0.063 | 0.070 | -0.903 | 11 | 0.386 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | pre_postpost | -0.056 | 0.058 | -0.967 | 11 | 0.354 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | pre_postpost | -0.003 | 0.004 | -0.741 | 11 | 0.474 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | pre_postpost | -0.006 | 0.025 | -0.238 | 11 | 0.816 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | pre_postpost | -0.041 | 0.047 | -0.885 | 11 | 0.395 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | pre_postpost | 0.631 | 0.511 | 1.234 | 11 | 0.243 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | pre_postpost | -4.914 | 4.239 | -1.159 | 11 | 0.271 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | pre_postpost | -0.076 | 0.433 | -0.175 | 11 | 0.864 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | pre_postpost | 2.985 | 0.809 | 3.690 | 11 | 0.004 | ** |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | pre_postpost | 2.506 | 0.771 | 3.250 | 11 | 0.008 | ** |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | pre_postpost | 0.163 | 0.154 | 1.055 | 11 | 0.314 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | pre_postpost | 0.061 | 0.028 | 2.152 | 11 | 0.054 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | pre_postpost | -0.018 | 0.018 | -1.004 | 11 | 0.337 | ns |
| x31_cd14_monocytes | fixed | NA | pre_postpost | -2.571 | 4.172 | -0.616 | 11 | 0.550 | ns |
| x32_cd16_non_classical_mono | fixed | NA | pre_postpost | -0.594 | 0.483 | -1.231 | 11 | 0.244 | ns |
| x33_cd16_classical_mono | fixed | NA | pre_postpost | -1.970 | 3.736 | -0.527 | 11 | 0.608 | ns |
| x34_hladr_cd56 | fixed | NA | pre_postpost | -2.146 | 3.763 | -0.570 | 11 | 0.580 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | pre_postpost | -0.078 | 0.060 | -1.301 | 11 | 0.220 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | pre_postpost | -1.310 | 3.138 | -0.417 | 11 | 0.684 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | pre_postpost | -1.292 | 1.590 | -0.812 | 11 | 0.434 | ns |
| x38_cd16_nk_cells | fixed | NA | pre_postpost | -0.601 | 0.521 | -1.153 | 11 | 0.273 | ns |
| x40_cd14_mdsc_mono | fixed | NA | pre_postpost | -1.314 | 2.014 | -0.652 | 11 | 0.528 | ns |
# extract statistically significant cells
cells_trtR_psig <- cells_trtR_coefonly %>%
filter(cells_trtR_coefonly$p.value < 0.05)
print(cells_trtR_psig$cell_type)
## [1] "x25_cd19_cd3_b_cells" "x26_cd27_ig_d_naive_b_cells"
sigcells_R <- cells_trtR_psig$cell_type
cells_R_long_edited %>%
filter(cell_type %in% sigcells_R) %>%
ggplot(aes(x = pre_post, y = cell_value, color = patient_id)) +
geom_line(aes(group = patient_id)) +
labs(x = "",
y = "",
title = "") +
theme_classic() +
facet_wrap(~ cell_type, scales = "free_y")
cells_R_long_edited %>%
filter(cell_type %in% sigcells_R) %>%
ggplot(aes(x = pre_post, y = cell_value, fill = pre_post)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(values = c("pre" = "gray",
"post" = "tomato")) +
geom_jitter() +
labs(x = "Pre vs post",
y = "Cell value",
title = "Significantly different cell types: pre vs post Red") +
facet_wrap(~ cell_type, ncol = 3, nrow = 2, scales = "free_y") +
theme_bw()
Naive B-cell pop (#26) significantly increase post-Red intervention, and also increase post-Yellow intervention. This could suggest there is a tomato effect.
I need to figure out how to make a function that extracts makes a QQplot for each model. But since IL-5 was significant, let’s look at how the residuals look for this cytokine (for every comparison) before log transforming
#qqplot_fx <- function(cytokine){
#qqnorm(resid(cyto_seq_models$cytokine))
#qqline(resid(cyto_seq_models$cytokine))
#}
#qqplot_fx("il_5")
#cyto_seq_results$
qqnorm(resid(cells_seq_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_seq_models$x25_cd19_cd3_b_cells))
Doesn’t look bad
qqnorm(resid(cells_trt_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trt_models$x25_cd19_cd3_b_cells))
doesn’t look as bad as sequence
qqnorm(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
qqnorm(resid(cells_trtY_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtY_models$x25_cd19_cd3_b_cells))
I will log transform the data and perform mixed linear modeling and paired t-tests for every comparison to see how the data compares.
Wrangle
cells_log2_long <- cells_long %>%
mutate(log2_cell_value = log2(cell_value)) %>%
filter(log2_cell_value != -Inf) # remove 0s
# lmer model function
cells_log2_seq_function <- function(cells_log2_long) {
lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_long, REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_seq_models <- map(split(cells_log2_long,
cells_log2_long$cell_type),
cells_log2_seq_function)
# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_seq_results <- map_df(cells_log2_seq_models,
tidy,
.id = "cell_type")
cells_log2_seq_results_coef <- cells_log2_seq_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
kable(cells_log2_seq_results_coef, format = "markdown", digits = 3)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | sequenceY_R | 0.002 | 0.035 | 0.064 | 10.000 | 0.950 | ns |
| x02_cd45_cd66b_grans | fixed | NA | sequenceY_R | 0.084 | 0.876 | 0.096 | 10.000 | 0.926 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | sequenceY_R | 0.226 | 0.317 | 0.712 | 10.000 | 0.492 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | sequenceY_R | 0.213 | 0.340 | 0.625 | 10.000 | 0.546 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | sequenceY_R | 0.573 | 0.367 | 1.562 | 10.000 | 0.149 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | sequenceY_R | 0.658 | 0.858 | 0.767 | 10.000 | 0.461 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | sequenceY_R | 0.500 | 0.516 | 0.969 | 10.000 | 0.355 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | sequenceY_R | 0.891 | 0.711 | 1.254 | 10.000 | 0.239 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | sequenceY_R | -0.186 | 0.431 | -0.432 | 10.000 | 0.675 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | sequenceY_R | 0.068 | 0.398 | 0.172 | 10.000 | 0.867 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | sequenceY_R | 0.116 | 0.507 | 0.230 | 10.000 | 0.823 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | sequenceY_R | 0.489 | 0.725 | 0.675 | 10.000 | 0.515 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | sequenceY_R | 0.055 | 0.515 | 0.108 | 10.000 | 0.916 | ns |
| x14_cd45ro_cd45ra | fixed | NA | sequenceY_R | 0.383 | 0.379 | 1.013 | 10.000 | 0.335 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | sequenceY_R | 0.964 | 0.886 | 1.088 | 10.000 | 0.302 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | sequenceY_R | 0.141 | 0.408 | 0.345 | 10.000 | 0.737 | ns |
| x17_cd25_cd127_tregs | fixed | NA | sequenceY_R | 0.334 | 0.449 | 0.745 | 10.000 | 0.473 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | sequenceY_R | 0.495 | 0.448 | 1.105 | 10.001 | 0.295 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | sequenceY_R | 0.727 | 0.688 | 1.056 | 10.000 | 0.316 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | sequenceY_R | 0.198 | 0.364 | 0.545 | 10.000 | 0.598 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | sequenceY_R | 0.422 | 0.442 | 0.956 | 10.000 | 0.362 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | sequenceY_R | 1.429 | 1.201 | 1.190 | 10.000 | 0.262 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | sequenceY_R | 1.660 | 1.442 | 1.151 | 10.000 | 0.276 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | sequenceY_R | 1.479 | 1.223 | 1.209 | 10.000 | 0.254 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | sequenceY_R | -0.074 | 0.456 | -0.162 | 10.000 | 0.874 | ns |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | sequenceY_R | 0.012 | 0.502 | 0.023 | 10.000 | 0.982 | ns |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | sequenceY_R | -0.461 | 0.551 | -0.837 | 10.000 | 0.422 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | sequenceY_R | -0.400 | 0.709 | -0.564 | 10.000 | 0.585 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | sequenceY_R | -0.795 | 0.652 | -1.220 | 9.743 | 0.251 | ns |
| x31_cd14_monocytes | fixed | NA | sequenceY_R | -0.405 | 0.448 | -0.905 | 10.000 | 0.387 | ns |
| x32_cd16_non_classical_mono | fixed | NA | sequenceY_R | -0.531 | 0.672 | -0.790 | 10.000 | 0.448 | ns |
| x33_cd16_classical_mono | fixed | NA | sequenceY_R | -0.353 | 0.447 | -0.789 | 10.000 | 0.449 | ns |
| x34_hladr_cd56 | fixed | NA | sequenceY_R | -0.354 | 0.442 | -0.800 | 10.000 | 0.442 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | sequenceY_R | -0.757 | 0.625 | -1.211 | 10.000 | 0.254 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | sequenceY_R | -0.278 | 0.457 | -0.608 | 10.000 | 0.557 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | sequenceY_R | -0.362 | 0.644 | -0.562 | 10.000 | 0.587 | ns |
| x38_cd16_nk_cells | fixed | NA | sequenceY_R | 0.213 | 0.631 | 0.338 | 10.000 | 0.742 | ns |
| x40_cd14_mdsc_mono | fixed | NA | sequenceY_R | -0.711 | 0.789 | -0.901 | 46.000 | 0.372 | ns |
| x41_cd66b_mdsc_grans | fixed | NA | sequenceY_R | -0.441 | 1.071 | -0.412 | 10.000 | 0.689 | ns |
# extract statistically significant cytokines
cells_log2_seq_psig <- cells_log2_seq_results_coef %>%
filter(cells_log2_seq_results_coef$p.value < 0.05)
print(cells_log2_seq_psig$cell_type)
## character(0)
*No sequence effects
# subset data into yellow results
cells_log2_Y_long <- cells_log2_long %>%
filter(intervention == "Yellow")
# lmer model function
cells_log2_Y_function <- function(cells_log2_Y_long) {
lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_Y_long, REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_Y_models <- map(split(cells_log2_Y_long,
cells_log2_Y_long$cell_type),
cells_log2_Y_function)
# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_Y_results <- map_df(cells_log2_Y_models,
tidy,
.id = "cell_type")
cells_log2_Y_results_coef <- cells_log2_Y_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
kable(cells_log2_Y_results_coef, format = "markdown", digits = 3)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | sequenceY_R | 0.053 | 0.049 | 1.086 | 10 | 0.303 | ns |
| x02_cd45_cd66b_grans | fixed | NA | sequenceY_R | -0.359 | 0.994 | -0.362 | 10 | 0.725 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | sequenceY_R | -0.072 | 0.320 | -0.226 | 10 | 0.826 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | sequenceY_R | -0.109 | 0.349 | -0.311 | 10 | 0.762 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | sequenceY_R | 0.488 | 0.392 | 1.245 | 10 | 0.242 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | sequenceY_R | 0.001 | 0.946 | 0.001 | 10 | 1.000 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | sequenceY_R | 0.160 | 0.533 | 0.299 | 10 | 0.771 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | sequenceY_R | 1.258 | 0.713 | 1.764 | 10 | 0.108 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | sequenceY_R | -0.085 | 0.563 | -0.151 | 10 | 0.883 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | sequenceY_R | 0.366 | 0.510 | 0.717 | 10 | 0.490 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | sequenceY_R | -0.325 | 0.518 | -0.627 | 10 | 0.544 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | sequenceY_R | -0.163 | 0.776 | -0.210 | 10 | 0.838 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | sequenceY_R | -0.347 | 0.528 | -0.658 | 10 | 0.525 | ns |
| x14_cd45ro_cd45ra | fixed | NA | sequenceY_R | 0.469 | 0.424 | 1.107 | 10 | 0.294 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | sequenceY_R | 1.360 | 0.860 | 1.582 | 10 | 0.145 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | sequenceY_R | 0.204 | 0.515 | 0.396 | 10 | 0.700 | ns |
| x17_cd25_cd127_tregs | fixed | NA | sequenceY_R | -0.016 | 0.504 | -0.033 | 10 | 0.975 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | sequenceY_R | 0.248 | 0.519 | 0.478 | 10 | 0.643 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | sequenceY_R | 0.540 | 0.754 | 0.716 | 10 | 0.490 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | sequenceY_R | 0.032 | 0.450 | 0.070 | 10 | 0.945 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | sequenceY_R | 0.140 | 0.497 | 0.281 | 10 | 0.784 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | sequenceY_R | 2.235 | 1.263 | 1.769 | 10 | 0.107 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | sequenceY_R | 1.188 | 1.581 | 0.752 | 10 | 0.470 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | sequenceY_R | 0.726 | 1.201 | 0.604 | 10 | 0.559 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | sequenceY_R | -0.028 | 0.425 | -0.066 | 10 | 0.949 | ns |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | sequenceY_R | 0.000 | 0.471 | 0.000 | 10 | 1.000 | ns |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | sequenceY_R | -0.228 | 0.570 | -0.400 | 10 | 0.698 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | sequenceY_R | -0.040 | 0.706 | -0.057 | 10 | 0.956 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | sequenceY_R | -0.650 | 0.737 | -0.881 | 10 | 0.399 | ns |
| x31_cd14_monocytes | fixed | NA | sequenceY_R | 0.164 | 0.489 | 0.334 | 10 | 0.745 | ns |
| x32_cd16_non_classical_mono | fixed | NA | sequenceY_R | 0.448 | 0.711 | 0.630 | 10 | 0.543 | ns |
| x33_cd16_classical_mono | fixed | NA | sequenceY_R | 0.217 | 0.498 | 0.436 | 10 | 0.672 | ns |
| x34_hladr_cd56 | fixed | NA | sequenceY_R | 0.144 | 0.482 | 0.298 | 10 | 0.772 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | sequenceY_R | -0.232 | 0.599 | -0.387 | 10 | 0.707 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | sequenceY_R | 0.258 | 0.507 | 0.509 | 10 | 0.622 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | sequenceY_R | -0.381 | 0.724 | -0.526 | 10 | 0.610 | ns |
| x38_cd16_nk_cells | fixed | NA | sequenceY_R | 0.454 | 0.679 | 0.670 | 10 | 0.518 | ns |
| x40_cd14_mdsc_mono | fixed | NA | sequenceY_R | -0.114 | 0.923 | -0.123 | 10 | 0.904 | ns |
# extract statistically significant cytokines
cells_log2_Y_psig <- cells_log2_Y_results_coef %>%
filter(cells_log2_Y_results_coef$p.value < 0.05)
print(cells_log2_Y_psig$cell_type)
## character(0)
Now there arent any significant cells for yellow after log transforming
# subset data into yellow results
cells_log2_R_long <- cells_log2_long %>%
filter(intervention == "Red")
# lmer model function
cells_log2_R_function <- function(cells_log2_R_long) {
lmer(log2_cell_value ~ sequence + (1|patient_id), data = cells_log2_R_long, REML = TRUE)
}
# break data up into subsets based on cytokine, then apply funtion to each subset
cells_log2_R_models <- map(split(cells_log2_R_long,
cells_log2_R_long$cell_type),
cells_log2_R_function)
# create single data frame of results. apply tidy from broom.mixed package to clean up results
cells_log2_R_results <- map_df(cells_log2_R_models,
tidy,
.id = "cell_type")
cells_log2_R_results_coef <- cells_log2_R_results %>%
filter(effect == "fixed") %>%
filter(term != "(Intercept)") %>%
add_significance(p.col = "p.value",
output.col = "p.sig",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
symbols = c("****", "***", "**", "*", "ns")
)
kable(cells_log2_R_results_coef, format = "markdown", digits = 3)
| cell_type | effect | group | term | estimate | std.error | statistic | df | p.value | p.sig |
|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | fixed | NA | sequenceY_R | -0.049 | 0.032 | -1.539 | 10.00 | 0.155 | ns |
| x02_cd45_cd66b_grans | fixed | NA | sequenceY_R | 0.527 | 0.808 | 0.652 | 10.00 | 0.529 | ns |
| x03_cd3_cd45_cd3_t_cells | fixed | NA | sequenceY_R | 0.524 | 0.345 | 1.519 | 10.00 | 0.160 | ns |
| x04_tc_rgd_cd3_ab_t_cells | fixed | NA | sequenceY_R | 0.534 | 0.365 | 1.463 | 10.00 | 0.174 | ns |
| x05_cd4_cd8_cd8_t_cells | fixed | NA | sequenceY_R | 0.658 | 0.367 | 1.792 | 10.00 | 0.103 | ns |
| x06_cd45ro_cd45ra_naive_cd8 | fixed | NA | sequenceY_R | 1.315 | 0.788 | 1.670 | 10.00 | 0.126 | ns |
| x07_cd46ro_cd45ra_cm_cd8 | fixed | NA | sequenceY_R | 0.841 | 0.543 | 1.550 | 10.00 | 0.152 | ns |
| x08_cd45ro_cd45ra_em_cd8 | fixed | NA | sequenceY_R | 0.525 | 0.743 | 0.707 | 10.00 | 0.496 | ns |
| x09_cd45r0_cd45ra_te_cd8 | fixed | NA | sequenceY_R | -0.287 | 0.318 | -0.903 | 10.00 | 0.388 | ns |
| x10_cd38_hladr_activated_cd8 | fixed | NA | sequenceY_R | -0.229 | 0.346 | -0.662 | 10.00 | 0.523 | ns |
| x11_cd4_cd8_cd4_t_cells | fixed | NA | sequenceY_R | 0.558 | 0.527 | 1.058 | 10.00 | 0.315 | ns |
| x12_cd45ro_cd45ra_naive_cd4 | fixed | NA | sequenceY_R | 1.141 | 0.723 | 1.578 | 10.00 | 0.146 | ns |
| x13_cd45ro_cd45ra_cm_cd4 | fixed | NA | sequenceY_R | 0.458 | 0.533 | 0.860 | 10.00 | 0.410 | ns |
| x14_cd45ro_cd45ra | fixed | NA | sequenceY_R | 0.298 | 0.369 | 0.806 | 10.00 | 0.439 | ns |
| x15_cd45ro_cd45ra_te_cd4 | fixed | NA | sequenceY_R | 0.568 | 1.041 | 0.546 | 10.00 | 0.597 | ns |
| x16_cd38_hladr_activated_cd4 | fixed | NA | sequenceY_R | 0.078 | 0.376 | 0.207 | 10.00 | 0.840 | ns |
| x17_cd25_cd127_tregs | fixed | NA | sequenceY_R | 0.685 | 0.418 | 1.638 | 10.00 | 0.132 | ns |
| x18_ccr4_cd4_total_ccr4_treg | fixed | NA | sequenceY_R | 0.743 | 0.410 | 1.813 | 10.00 | 0.100 | ns |
| x19_cd45ra_cd45ro_ccr4_treg_naive | fixed | NA | sequenceY_R | 0.914 | 0.860 | 1.062 | 10.00 | 0.313 | ns |
| x20_hladr_total_ccr4_treg_activated | fixed | NA | sequenceY_R | 0.365 | 0.324 | 1.126 | 10.00 | 0.286 | ns |
| x21_cd45ra_cd45ro_ccr4_treg_memory | fixed | NA | sequenceY_R | 0.705 | 0.409 | 1.724 | 10.00 | 0.115 | ns |
| x22_cxcr3_ccr6_th1 | fixed | NA | sequenceY_R | 0.623 | 1.191 | 0.523 | 10.00 | 0.613 | ns |
| x23_cxcr3_ccr6_th2 | fixed | NA | sequenceY_R | 2.133 | 1.316 | 1.620 | 10.00 | 0.136 | ns |
| x24_cxcr3_ccr6_th17 | fixed | NA | sequenceY_R | 2.233 | 1.285 | 1.738 | 10.00 | 0.113 | ns |
| x25_cd19_cd3_b_cells | fixed | NA | sequenceY_R | -0.120 | 0.506 | -0.237 | 10.00 | 0.818 | ns |
| x26_cd27_ig_d_naive_b_cells | fixed | NA | sequenceY_R | 0.023 | 0.565 | 0.041 | 10.00 | 0.968 | ns |
| x27_cd27_ig_d_memory_b_cells | fixed | NA | sequenceY_R | -0.695 | 0.554 | -1.253 | 10.00 | 0.239 | ns |
| x28_cd27_ig_d_memory_resting_b_cells | fixed | NA | sequenceY_R | -0.759 | 0.741 | -1.025 | 10.00 | 0.330 | ns |
| x30_cd27_cd38_plasmablasts | fixed | NA | sequenceY_R | -0.797 | 0.531 | -1.503 | 8.25 | 0.170 | ns |
| x31_cd14_monocytes | fixed | NA | sequenceY_R | -0.974 | 0.465 | -2.096 | 10.00 | 0.063 | ns |
| x32_cd16_non_classical_mono | fixed | NA | sequenceY_R | -1.510 | 0.733 | -2.061 | 10.00 | 0.066 | ns |
| x33_cd16_classical_mono | fixed | NA | sequenceY_R | -0.922 | 0.458 | -2.013 | 10.00 | 0.072 | ns |
| x34_hladr_cd56 | fixed | NA | sequenceY_R | -0.851 | 0.459 | -1.854 | 10.00 | 0.093 | ns |
| x35_cd16_cd123_cd11c_p_dc | fixed | NA | sequenceY_R | -1.282 | 0.684 | -1.874 | 10.00 | 0.090 | ns |
| x36_cd16_cd123_cd11c_m_dc | fixed | NA | sequenceY_R | -0.813 | 0.463 | -1.757 | 10.00 | 0.110 | ns |
| x37_cd56_cd161_cd123_nk_cells | fixed | NA | sequenceY_R | -0.343 | 0.671 | -0.512 | 10.00 | 0.620 | ns |
| x38_cd16_nk_cells | fixed | NA | sequenceY_R | -0.028 | 0.629 | -0.044 | 10.00 | 0.965 | ns |
| x40_cd14_mdsc_mono | fixed | NA | sequenceY_R | -1.309 | 0.717 | -1.825 | 10.00 | 0.098 | ns |
| x41_cd66b_mdsc_grans | fixed | NA | sequenceY_R | -0.441 | 1.071 | -0.412 | 10.00 | 0.689 | ns |
# extract statistically significant cytokines
cells_log2_R_psig <- cells_log2_R_results_coef %>%
filter(cells_log2_R_results_coef$p.value < 0.05)
print(cells_log2_R_psig$cell_type)
## character(0)
Now there arent any significant cells for red after log transforming
# before log transform
qqnorm(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_trtR_models$x25_cd19_cd3_b_cells))
# after log transform
qqnorm(resid(cells_log2_R_models$x25_cd19_cd3_b_cells))
qqline(resid(cells_log2_R_models$x25_cd19_cd3_b_cells))
Looks even worse after log transform